n_tasks = 9;
pca_dimensionality = 30;
n_training = 7;
n_sampleruns = 2;

load('faces.mat');

% compute mean vector and PCA projection from all data
x = [];
for img=1:n_tasks;
  x=[x;males{img}];
end;
for img=1:n_tasks;
  x=[x;females{img}];
end;
datamean=mean(x,1);
datacov=cov(x);
[pcaproj,D]=eig(datacov);
d=diag(D); [s,I]=sort(-d);
pcaproj=pcaproj(:,I(1:pca_dimensionality));

%[pcaproj,data] = princomp(x);
%pcaproj=pcaproj(:,1:pca_dimensionality);


% transform data into format required by RSL code: [features 1 classlabel] for each image
tasksmales = cell(n_tasks,1);
tasksfemales = cell(n_tasks,1);
for i=1:n_tasks,
  % Substract overall data mean,
  % apply PCA,
  % append a constant feature (value 1) to the end of the data,
  % and append the class (here -1 males, +1 females)
  tempdata1=males{i};
  tempdata1=[(tempdata1-repmat(datamean,[size(tempdata1,1) 1]))*pcaproj ones(size(tempdata1,1),1) -1*ones(size(tempdata1,1),1)];
  tasksmales{i} = tempdata1;
  tempdata2=females{i};
  tempdata2=[(tempdata2-repmat(datamean,[size(tempdata1,1) 1]))*pcaproj ones(size(tempdata2,1),1) 1*ones(size(tempdata2,1),1)];
  tasksfemales{i} = tempdata2;
end;

% create a stratified random permutation of the 40 persons;
% leave 5 males and 5 females in training sets, rest in test sets
outputfilter=cell(3,n_tasks,n_training,n_sampleruns);
avg_error = cell(3,n_tasks,n_training);
for i=1:numel(avg_error)
    avg_error{i} = 0;
end


% ----------
% Main Loop
% ----------
stlqualities = zeros(n_training,n_tasks);
poolqualities = zeros(n_training,n_tasks);
rslqualities = zeros(n_training,n_tasks);
total_stlqualities = zeros(n_training,n_tasks);
total_poolqualities = zeros(n_training,n_tasks);
total_rslqualities = zeros(n_training,n_tasks);
for n_run=1:n_sampleruns
    for n_trainpersons=1:n_training
        prev_error = avg_error;
        traintasks=cell(n_tasks,1);
        testtasks=cell(n_tasks,1);
        for i=1:n_tasks,
            ord1 = randperm(20);
            ord2 = randperm(20);
            traintasks{i}=[tasksmales{i}(ord1(1:n_trainpersons),:); tasksfemales{i}(ord2(1:n_trainpersons),:)];
            testtasks{i}=[tasksmales{i}(ord1((n_trainpersons+1):end),:); tasksfemales{i}(ord2((n_trainpersons+1):end),:)];
        end;
        
        % STL
        temp_acc = 0;
        for task=1:n_tasks,
            [outputfilter{1,task,n_trainpersons,n_run} stlqualities(n_trainpersons,task) ] = run_vb_stl( traintasks{task},testtasks{task},pca_dimensionality );
            
        end;
       
      
        % pooled STL
        for task=1:n_tasks,
          % create a pooled data set
            tempdata = [];
            for t=1:n_tasks,
                tempdata=[tempdata;traintasks{t}];
            end;
            [outputfilter{2,task,n_trainpersons,n_run}   poolqualities(n_trainpersons,task)] = run_vb_pooled( tempdata,testtasks{task},pca_dimensionality );
            
        end;
       
     
        
        % RSL
        for task=1:n_tasks,
            [outputfilter{3,task,n_trainpersons,n_run } rslqualities(n_trainpersons,task)] = run_vb_rsl(traintasks,testtasks{task},task,pca_dimensionality);
           
        end;
       
        
    end
    total_rslqualities = total_rslqualities + rslqualities;
    total_poolqualities = total_poolqualities + poolqualities;
    total_stlqualities = total_stlqualities + stlqualities;
    % Add previous error
    %avg_error = CellOperation.add(avg_error,prev_error);
end

% calculate average error
%avg_error = CellOperation.div(avg_error,n_sampleruns);

total_rslqualities = total_rslqualities ./ n_sampleruns;
total_poolqualities = total_poolqualities ./ n_sampleruns;
total_stlqualities = total_stlqualities ./ n_sampleruns;


total_rslqualities
total_poolqualities
total_stlqualities


figure;
for method=1:3,
  for task=1:n_tasks,
    subplot(3,n_tasks,(method-1)*n_tasks+task);
    imagesc(reshape(pcaproj*outputfilter{method,task,1,1},[16 16]));
    colormap(gray);
  end;
end;

figure;
for method=1:3,
  for task=1:n_tasks,
    subplot(3,n_tasks,(method-1)*n_tasks+task);
    imagesc(reshape(pcaproj*outputfilter{method,task,7,1},[16 16]));
    colormap(gray);
  end;
end;

